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Abstract. An abstract should be given 

Semi-empirical molecular dynamics is used to simulate several gaseous atomic hy- 
drogen interactions with hydrocarbon grains in space: recoil, adsorption, diffusion, 
chemisorption and recombination into molecular hydrogen. Their probabilities are deter- 



q . mined as a function of initial velocity of gas atoms. The equilibrium hydrogen coverage 
H ' 

of free grains is then derived. These data can be used in calculations of material and 
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energy balance as well as rates of hydrogen recombination on grains. The role of grain 
temperature is also considered. 
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1. Introduction 

The primary purpose of the present work is to study several processes which occur when 
H atoms impinge on carbonaceous grains, and thus bear upon the recombination of these 
atoms into diatomic molecules. This, in turn, affects chemical and energetic equilibria in 
gas and grains. Here, focus is set on such equilibria in environments typical of carbon- 
rich PDRs (Photon-Dominated Regions). Gas-surface interaction has been the subject of 
intensive and extensive experimental work (e.g. Winkler and Rendulik 1992, Donnet et 
al. lT§9"5|l as well as an abundant literature, including several treatises (see Bortolani et al. 
119941 Masel IT99f)[) . In general, however, available information relates to very low target 
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temperatures (e.g. Pirronello et al. [1999 Katz et al. 1 19991 Cazeaux and Tielens 2004) 
and/or pure and regular dielectric or metallic surfaces of technological interest (see many 
issues of Surface Science). Here, instead, we are interested in hydrogenated, amorphous, 
semiconducting carbon surfaces. 

When an H atom hits such a surface, it may recoil and be scattered away. Or it 
may be attracted and captured by a dangling C-bond (dehydrogenated site) and stick 
to the surface (strong chemisorption) . Otherwise, it may interact electronically with one 
particular H atom of the surface, so strongly as to form an H2 molecule, which promptly 
leaves the surface. This process is called abstraction or E-R (Elcy-Rideal) process; it can 
be direct or indirect (retarded) according to whether it takes place promptly or after 
some wandering over the surface (see Stecher and Williams 1966). Such surface diffusion 
is induced by the initial velocity of the projectile and enhanced by thermal agitation of 
the target atoms. 

If the H projectile has a low initial velocity, it may be adsorbed in a shallow (ph- 
ysisorption) well, making a weak bond to the surface, which allows diffusion. When the 
grain temperature is not too high, the atom may stay on the surface for a long time. 
Then, a second projectile may follow a similar course and diffuse to meet the first one, 
in which case a molecule is formed and leaves the surface: this is the H-L (Hinshelwood- 
Langmuir) process. Both the E-R and H-L mechanisms contribute to the recombination 
rate R (cm 3 s _1 ; see Duley and Williams 11984(1 . in various amounts, according to gas 
density and grain temperature. 

Finally, if the projectile velocity is very high, it may kick out the target H atom 
{expulsion) and, sometimes, substitute to it at the same C site {substitution). The hy- 
drogen coverage of the grain surface, and hence the recombination rate as well as the IR 
(InfraRed) spectrum, strongly depend on the balance of all these processes, which may 
all occur on the same grain when the gas atom velocities cover a wide range. 

All these processes are sensitive to velocity, angle and site of impact; hence the great 
difficulty of the experimental approach. As to the classical theoretical approach, it re- 
quires tailoring essential parameters (e.g. activation energies; see Duley 1199 6. Habart 
et al. |2004) to fit observations. Parneix and Brechignac (1998), instead, used numerical 
simulation to study H2 recombination on graphite from first principles. Here, Molecular 
Dynamics is also used to simulate numerically not one particular, but all the above 
phenomena occuring simultaneously. For this purpose, the target is not chosen to be a 
surface, but a finite, small grain, and not one, but several projectiles are launched at the 
same time towards it from different directions, with the same velocity, Vo, set, at the 
start of each run, at a different value, between 1 and 500 km.s_i (which includes the 
range of gas temperatures of interest, between and 1000 K). Also, because the grain 
temperature is much higher in PDRs than in molecular clouds, the model target was 
given initial temperatures ranging from to 365 K. One essential feature of this simu- 
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lation is the fact that the target structure is not rigid: both the nuclear and electronic 
positions change under the electrical influence of the incoming H atom, allowing energy 
exchange between them; otherwise, only elastic recoil is possible. More generally, all the 
physics involved is included a priory in the chemical code which governs the dynamics, 
and no assumption is made on any of the parameters of interest (e.g. sticking coefficient 
on, or H coverage of, grains, or H2 photo-dissociation rate or activation energies). 

2. The chemical code 

A description of the software package used for simulation can be found in Papoular ( 2001 ) . 
The particular code used here is AMI (Austin Model 1), as proposed, and optimized for 
carbon-rich molecules, by Dewar et al. (1985). This semi-empirical method combines a 
rigorous quantum-mechanical formalism with empirical parameters obtained from com- 
parison with experimental results. It computes approximate solutions of Schroedinger's 
equation, and uses experimental data only when the Q.M. calculations are too difficult 
or too lengthy. This makes it more accurate than poor ab initio methods, and faster than 
any of the latter. 

The molecular dynamics relies on the Born-Oppenheimer approximation to determine 
the motions of atoms under nuclear and electronic forces due to their environment. At 
every step, all the system parameters are memorized as snapshots so that, after comple- 
tion of the run, a movie of the reactions can be viewed on the screen, and the trajectory 
of any atom followed all along. This makes for a better understanding of the details of 
mechanisms and outcomes. Note that the dynamics of bond dissociations and formation 
can only be simulated by using Unrestricted Hartree-Fock (UHF) wave functions in the 
Q.M. part of the calculation (see Szabo and Ostlund ll989f) . The elementary calculation 
step was set at 10~ 15 s for initial projectile velocities smaller than 200 km.s -1 and 10~ 16 
above. 

Semi-empirical simulation methods do not account for purely quantum-mechanical 
phenomena like tunneling or zero-point energy. In the present context of relatively high 
grain temperature, the former is negligible as regards surface diffusion, and the latter 
hardly affects the probability and energetics of the physico-chemical reactions involved. 

3. The model dust grain 

This study is restricted to carbonaceous IS (InterStellar) dust. IR spectra of such dust 
reveals the presence of both aromatic and aliphatic C-C and C-H bonds (e.g. Pendleton 
et al. I1994f) . Graphitic dust has been abundantly modeled (see Parneix and Brechignac 
1998, Ko et al. 2002). Here, the focus is on amorphous aliphatic dust. In order to provide 
a model of fully hydrogenated surface (Fig. 1, right), one 3-D structure of 100 atoms was 
assembled from 66 C atoms and 34 H atoms, forming two 6-membered, one 5-membered 
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and one 4-membered rings, and several aliphatic chains. 11 H atoms are involved in CH 
sites or functionalities (tertiary carbon), 44 in CH2 sites (secondary carbon) and 9 in 
CH 3 sites (primary carbon) . The bond lengths and angles of this initial structure were 
then automatically optimized by the software, searching for the minimum overall poten- 
tial energy. The final system is enclosable in a box of size 15*12*7 AA 3 (the cluster 
sizes adopted here are a compromise between computational time and relevance to real 
IS dust). In a preliminary run, the system is brought up to the desired temperature, T, 
meaning that the atoms are given random velocities, increasing progressively in a prede- 
termined number of calculation steps, subject to bond constraints, up to a Maxwellian 
distribution characteristic of T. Such a set-up provides information on all the processes 
described above, except recombination of H atoms on dangling C-bonds, which are absent 
here. The latter is better studied by means of the second, completely de-hydrogenated, 
target, Fig. 1, left. It includes 33 C atoms, of which 2 C-sites with 3 dangling bonds, 10 
sites with 2, 11 with only 1 and 10 with no such bond. Its size is 9*5*4 A 3 . 

In order to study the energy balance, it is necessary to know the net result of each 
elementary process, which requires knowledge of the C-H and H-H bonding energies. 
These are obtained from the cold, optimized targets by extracting a particular H atom 
and noting the increase in overall potential energy. They are thus found to be 109.4 
kcal/mole for H2, 107 for aromatic C-H and about 96 for all aliphatic C-H bonds. The 
corresponding values in eV are 4.75, 4.65 and 4.15 respectively, and the bond lengths, 
0.68, 1.1 and about 1.1 AA , respectively. By contrast, a weak chemical bond, 2.1 eV, is 
formed with one C atom of a graphitic plane by an H atom, 1.15 A above the C atom, 
which is itself slightly uplifted in the process. 

Physisorption sites have been found by minimizing the potential energy of an H atom 
deposited at some arbitrary, but typical, locations around and at a few angstroms from 
the target. For the hydrogenated aliphatic grain, most of them are of 2 types: a) at 
2.85 A from a C atom, along the C-H bond, the potential well depth or bond energy 
is 45 meV, same as for an aromatic C-H site studied independently; b) between two 
adjacent H atoms, at about 2 A from both, the well depth is about 80 meV. For the 
dehydrogenated grain, physisorption sites are absent due to the long range attraction of 
the dangling bonds. 

In dynamic calculations, no activation barrier for abstraction was found in the exact 
direction of aliphatic C-H bonds, but there is an increasing potential barrier as the angle 
of incidence increases, and no abstraction was observed beyond 60 degrees. By contrast 
with aliphatic bonds and weak bonds above a graphitic plane, no abstraction at all was 
observed in the case of strong aromatic C-H bonds. This must be traced to the greater 
rigidity of the latter. Indeed, freezing the C-H bond (i.e. letting it retain the same length 
and direction during the calculation) suppresses abstraction even in aliphatic sites. The 
importance of bond elasticity is highlighted by the direct observation of the harpoon effect 
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Fig. 1. Right: the hydrogenated target with saturated C-bonds; 64 C atoms (larger cir- 
cles) and 36 H atoms (smaller circles); enclosing box 15*12*7 A 3 . Left: the pure-carbon 
target with dangling (free) bonds; 33 C atoms; 9*5*4 A 3 . The isolated H atoms sur- 
rounding the targets are the projectiles, 10 and 6 respectively, which are all given the 
same initial velocity, directed to a central C atom in each target. 

(see Atkinfl998) just before abstraction occurs: as the projectile approaches a target H 
atom, both their electron distributions change so that the attraction increases between 
them while the adjacent C-H bond expands. These changes are reversible until the saddle 
point, or transition state, is reached, when the distance between the H atoms is 1.3 A . 
Then, this distance is suddenly reduced by one half and the H2 molecule is formed and 
escapes. The saddle point clearly shows up in the potential curve, in the form of a small 
bump within the main well. No such phenomenon occurs with the aromatic C-H bond. 

4. Data analysis and results 

In a typical experiment, the target is surrounded by 10 or 6 H atoms (for the hydrogenated 
and dehydrogenated targets, respectively) at a distance larger than 4 A and as far as 
possible from each other, to minimize interference between them and occurences of L-H 
events between adsorbed projectiles. Each is given a velocity Vq, directed to an arbitrary, 
but central, atom of the target; Vo spans the range between 1 and 500 A ps" 1 (10 4 to 
5 10 6 cm.s -1 ) corresponding to impact energies from 5 10 -5 to 12.5 eV. The calculation 
is then started and the trajectory of each projectile and target atom is followed in time 
and visualized on the computer screen. The calculation extends as long as necessary for 
all projectiles to interact (almost) completely: this usually takes less than 2 ps (2 10~ 12 
s). 

Let N r be the number of completed reactions of type r between a projectile and an 
H or C site of which s r are available out of a total of s target sites. A corresponding 
coverage is defined as c r — s r /s. If a single projectile were used and the experiment 
repeated a large number of times, p, so that N r were very large, a probability could be 
defined for reaction r by the first of the two equalities 
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Fig. 2. Figure 2a: recombination probability coefficent,fc c , for the pure -C target at 
K initial temperature. Figure 2b: probability coefficients for the H-saturated target at 
K: prompt or direct abstraction (fc a , heavy full line); transient physisorption or indirect 
abstraction (k p , thin); expulsion (fc e , dots). Figure 2c: Hydrogen coverage, Ch, of a grain 
in equilibrium at K (dots) and 350 K (crosses). Figure 2d: the recombination rate, R, 
on a grain in equilibrium at K (dots) and 350 K (crosses). All fc r 's, and hence coverage 
and R, are functions of the initial velocity, Vo, assumed the same for all projectiles. Note 
the different scale for fig. 2d. 



p r = N r /p = k r .c r ; 

the second equality defines k r as a coefficient specific of this reaction which is, in 
general, a function of Vo and temperature T. Here, computer time is saved by using p 
projectiles simultaneously in each run instead of repeating the run p times and, in an 
attempt to capture the main features of the behaviour of k r for the various possible reac- 
tions, the crude assumption is made, that the notion of a probability p r may be extended 
to the present calculations. Now, that probability is, by definition, proportional to the 
number of relevant sites and to a corresponding effective crossection; the proportionality 
constant is found by rewriting the coverage as (s r /S).(S/s), where S is the surface area 
of the grain. Whence the effective crossection per such site, a r = k r .S/s, which gives a 
geometrical meaning to k r . Since the distance between two adjacent target sites is in the 
range 1 to 1.6 A , S/s may be estimated at 2 A 2 . 



R. Papoular: Gas-surface interactions 7 

Figure 2a represents the probability coefficient k c for recombination on the dchydro- 
genated target (at zero initial temperature), as deduced from the second of the equalities 
above, with c r =l and p=6. There seems to be a threshold for recombination at about 
5 10 4 cm.s -1 , indicating an activation energy of ^35 meV. This, however, is not associated 
with any bump in the static potential curve. It is due to the fact that the approaching 
projectile continuously gives away some of its kinectic energy to the deformable target. If 
and when it comes to a standstill, at ~ 4 A from the target, the latter relaxes and repels 
the projectile. Beyond velocity threshold, the latter retains all along enough energy to 
reach the point of no return, at the edge of the chemisorption well. 

Far beyond threshold, the argument is reversed: recombination is only possible if the 
projectile can give away enough of its energy to the target. As Vq increases, this becomes 
more and more difficult because the amount of energy to be evacuated increases while 
the time available to do that shortens. Hence, the recoil probability increases and recom- 
bination becomes rarer, while the minimum distance of approach before recoil decreases 
and tends to ~ 0.5 A . At intermediate velocities, optimal conditions for recombination 
are met and its probability reaches 1. The target temperature hardly affects these results: 
at 350 K, the velocity threshold for recombination is very slightly increased. 

Figure 2b represents the probability coefficients for (transient) physisorption on, and 
prompt (direct) abstraction and expulsion of H atoms from, the saturated target at zero 
initial temperature. Here, p=10 and c r is 1, 0.64 and 0.64 respectively for the 3 types of 
reaction. The threshold for expulsion is 250 A /ps and its probability increases quickly 
beyond. The corresponding threshold energy is ~ 3.1 eV, as compared with 4 eV bond 
energy; the difference is contributed by target deformation. 

At velocities lower than 20 A /ps, contrary to the dehydrogenated projectile 
can now be caught in the physisorption belt, and diffuse all around the target. During 
this journey, it ultimately abstracts an H atom from a hydrogenated site (delayed or 
indirect abstraction) or recombines on a dangling C bond, if there is any. Most often, 
this occurs in less than 2 ps because diffusion is very fast (>1 A /ps). At higher grain 
temperatures, the projectiles are more easily repelled before hitting the surface: at 350 K, 
the probability of physisorption is nearly uniformly reduced by a factor 2; also, diffusion 
is accelerated, but this temperature is not high enough for the adsorbate to be desorbed. 
In the following, we shall therefore assume that, for gas densities and temperatures of 
interest, recombination or indirect abstraction happens before another projectile stands 
a chance of being physisorbed and participating in a H-L recombination process. 

With these coefficients at hand, it is now possible to derive the hydrogen coverage of a 
free grain. Assuming exact balance between abstraction {k a with a delayed contribution 
by diffusing adsorbed atoms, k p ) and expulsion (fc e ) from hydrogenated sites, on the one 
hand, and recombination (k c with a delayed contribution by diffusing adsorbed atoms) 
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on the other hand, the equilibrium coverage is found to be ch = 1 — cc = k + fc C ^fc P +fc ■ 
Figure 2c represents this coverage for and 350 K initial temperatures. They differ no- 
tably only below 5 A /ps, because of the different probabilities of adsorption. By definition 
(see above) , knowledge of ch and k r finally allows one to determine the probability, p r , of 
each of the surface phenomena of interest, which, in turn, is necessary to determine the 
energy balance. Thus, the probability of recombination per H atom impact on a typical 
grain is: P — (k c + k p )cn- Now, the recombination rate is R = 1/3 Vb-PS, where S is 
the (carbonaceous) grain surface per H atom in space and is taken to be 2 10 -21 cm 2 . 
This quantity is represented in Fig. 2d, as a function of an effective gas temperature 
T e ff = ttitihVq /8fcs. In effect, the fc r 's should be replaced by integrals over the curves 
of Fig. 2a and 2b weighted by a distribution of velocities. Figures 2c and 2d show that 
grain temperature in this range has little relevance to H coverage or R. On the other 
hand, R is strongly dependent on gas temperature and dust abundance. Considering the 
large uncrtainties on the latter, the trend and range of R are in rough agreement with 
those deduced from observations by Habart et al. {2004 ). 

Regarding energy balance, note that, when recombination of an H atom occurs on a 
target dangling bond, all the available chemical energy (the C-H bond energy) is deposited 
in the grain, as is the initial projectile energy. By contrast, when abstraction occurs, the 
available energy is the sum of the projectile kinetic energy and the H-H bond energy but 
most of it goes into breaking the target C-H bond and feeding the internal and external 
kinetic energies of the newly formed H2 molecule. The detailed break-up depends on the 
particulars of the impact, but little kinetic energy is left in the grain and the molecule 
leaves promptly with a velocity higher than 70 A /ps. Similarly, in recoil, only a small 
fraction of the projectile kinetic energy is deposited in the target. 
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